\documentclass{report}

\usepackage{amsthm}
\usepackage{ulem}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{graphicx}

\usepackage{titlesec}
\usepackage{listings}
\usepackage{color}
\usepackage{wrapfig}
\usepackage{float}
\definecolor{dkgreen}{rgb}{0,0.6,0}
\definecolor{gray}{rgb}{0.5,0.5,0.5}
\definecolor{mauve}{rgb}{0.58,0,0.82}

\usepackage{url}

\titleformat{\chapter}[hang]
  {\normalfont\huge\bfseries}
  {\thechapter\ \ }
  {0pt}{}
  
\titlespacing*{\chapter}{0pt}{40pt}{20pt}
\titleformat*{\section}{\Large\bfseries}
\titleformat*{\subsection}{\large\bfseries}

\usepackage{etoolbox}
\makeatletter
\patchcmd{\chapter}{\if@openright\cleardoublepage\else\clearpage\fi}{\par}{}{}
\makeatother

\lstset{frame=tb,
  language=Java,
  aboveskip=3mm,
  belowskip=3mm,
  showstringspaces=false,
  columns=flexible,
  basicstyle={\small\ttfamily},
  numbers=none,
  numberstyle=\tiny\color{gray},
  keywordstyle=\color{blue},
  commentstyle=\color{dkgreen},
  stringstyle=\color{mauve},
  breaklines=true,
  breakatwhitespace=true
  tabsize=3
}

\usepackage{titling}
\newcommand{\subtitle}[1]{%
  \posttitle{%
    \par\end{center}
    \begin{center}\large#1\end{center}
    \vskip0.5em}%
}

\def\myprod#1#2{\prod_{\substack{#1\\#2}}}

\def\myvec#1#2{\begin{array}{c}#1\\#2\end{array}}

\def\mych#1#2{\left[ \myvec{#1}{#2} \right]}

\def\myth#1#2{\theta\left[ #1 \right] (#2,\tau)}

\def\mythp#1#2#3{\theta\left[ #1 \right]#2 (#3, \tau)}

\def\mythv#1#2#3{\myth{\myvec{#1}{#2}}{#3}}

\begin{document}

\nocite{*}

\title{Properties of Division Polynomials and Their Resultants}

\subtitle{Presentation for Masters Talk}

\date{\today}

\author{James D. Maxwell\\Department of Mathematics\\University of Colorado at Boulder}

\maketitle

\begin{abstract}
Let $E$ be an elliptic curve over $\mathbb{C}$, defined by $y^2 = x^3 + ax + b$.  
Let $\psi_n(x)$ be the associated division polynomial for an odd integer $n$. 
Let $\wp(z)$ be the Weierstrass $\wp$-function attached to $E$.
Define the derivation, $D$, on the function field of $E$ by
$D = 2y\frac{d}{dx}$.  Using this operator, we will computationally verify a formula
of Grant, that  
\[
res(\psi_n(x), D^2\psi_n(x))^2 = (-1)^{\frac{n-3}{2}}
2^{n^4-1}\pi^{n^4-n^2}n^{3n^2-3}\eta(\tau)^{2n^4-3n^2+1}Derivprod[n](\tau),
\]
where
\[ Derivprod[n](\tau) = \myprod{0 \leq a,b < n}{(a,b) \neq (0,0)}\mythp{1/2 +
\myvec{a/n}{b/n}}{'}{0},\]
$\theta$ is the Riemann theta function with
characteristic, and $res(\psi_n(x), D^2\psi_n(x))$ is the resultant of
$\psi_n(x)$ and $D^2\psi_n(x)$ with respect to $x$.

We will also similarly verify for various $n$ and $m$ that $res(\psi_n(x),
\psi_m(x))$ is a power of the discriminant $\Delta = -16(4a^3 + 27b^2)$, when $n$ and $m$ are relatively prime and odd, and that
$res(\frac{\psi_n(x)}{2y}, \psi_m(x))$ is a power of the discriminant when $n$
and $m$ are relatively prime, $n$ is even, and $m$ is odd.
\end{abstract}

\chapter{Introduction}

The intention of this paper is to calculate, using software I have written for this project, 
sample values for formulae related to elliptic curves, division polynomials, and
theta functions, 
including formulae from a paper written by David Grant (see \cite{Grant}).  
Grant's paper derives a result related to resultants and 
division polynomials of elliptic curves.  My paper is organized as follows.

Chapters one though five give a quick overview of elliptic curves - how they are 
defined, and their properties - division polynomials, resultants, and 
theta and Weierstrass functions related to an elliptic curve.

Chapter six is a summary of Grant's paper.  Chapter seven discusses my software
for elliptic curves, division polynomials, and related theta and
Weierstrass functions.
Chapter eight will verify some well known formulae related to elliptic
curves and Jacobi elliptic functions, and then conclude by testing Grant's
derivation by calculating some explicit values. 

\chapter{Elliptic Curves}

\section{Elliptic Curves Via Lattices}

We begin our discussion with a definition.

\newtheorem{mydef}{Definition}
\newtheorem{mythm}{Theorem}

\begin{mydef}
Let $\omega_1$, $\omega_2 \in \mathbb{C}$, be such that $\omega_1$ and
$\omega_2$ are $\mathbb{R}$-linearly independent.  We call $\Lambda = \mathbb{Z}\omega_1 \oplus \mathbb{Z}\omega_2$
a \emph{lattice} in $\mathbb{C}$. 
\end{mydef}

\begin{mydef}
A function, $f$, is \emph{doubly-periodic} if there is a lattice $\Lambda
\subseteq \mathbb{C}$ such that $f(z + \lambda) = f(z)$ for all $\lambda \in
\Lambda$.
\end{mydef}

\begin{mydef}
A function over $\mathbb{C}$ that is doubly-periodic and meromorphic on
$\mathbb{C}$ is referred to as an \emph{elliptic function}.
\end{mydef}

\begin{mydef}
Given a lattice, $\Lambda \subseteq \mathbb{C}$, we define the \emph{Weierstrass
$\wp$-function} attached to $\Lambda$ by \[ \wp(z, \Lambda) = \wp(z) =
\frac{1}{z^2} + \sum_{\substack{\omega \in \Lambda\\\omega \ne 0}} \left( \frac{1}{(z-\omega)^2} - \frac{1}{\omega^2} \right) \]
\end{mydef}

\begin{mydef}
Given a lattice, $\Lambda \subseteq \mathbb{C}$, define the \emph{Eisenstein
series} as \[ G_k(\Lambda) = G_k = \sum_{\substack{\omega \in \Lambda\\\omega \ne 0}}
\omega^{-k} \] where $k \geq 4$ is an even integer.
\end{mydef}

It is customary to assign $g_2 = 60G_4$ and $g_3 = 140G_6$.

\begin{mythm}
Given a lattice, $\Lambda \subseteq \mathbb{C}$, the Weierstrass $\wp$-function
attached to $\Lambda$ satisfies the differential equation
\[ \wp^\prime(z)^2 = 4\wp^3 - g_2\wp - g_3 \] 
\end{mythm}

\begin{proof}
See \cite{Wash}.
\end{proof}

\section{Elliptic Curves Over $\mathbb{C}$}
\begin{mydef}
An elliptic curve is defined by an equation of the form 
$y^2 = x^3 + ax + b$, where $a,b \in \mathbb{C}$ and $\Delta = -16\left(
4a^3 + 27b^2\right) \neq 0$.
\end{mydef}

So given a lattice,
$\Lambda = \mathbb{Z}\omega_1 \oplus \mathbb{Z}\omega_2$ in $\mathbb{C}$, 
there is a corresponding elliptic curve $E_{\Lambda}$ defined by $y^2 = x^3 -
\frac{g_2}{4}x - \frac{g_3}{4}$, and a surjective corresponding map
$(\mathbb{C}-\Lambda)/\Lambda \to E(\mathbb{C})$ defined by 
\begin{align}
 y &= \frac{\wp'(z)}{2}\\
 x &= \wp(z)
\end{align}

Conversely, given an elliptic curve defined by an equation $y^2 = x^3 + ax + b$, there exists
a lattice $\Lambda = \mathbb{Z}\omega_1 \oplus \mathbb{Z}\omega_2$ such that 
\begin{align}
a &= -\frac{g_2(\omega_1,\omega_2)}{4},\\
b &= -\frac{g_3(\omega_1,\omega_2)}{4}.
\end{align} 

This way we can define an elliptic curve from a lattice, or equivalently from an
equation of the form $y^2 = x^3 + ax + b$.

\section{Calculation of Periods}

Let $f = y^2 = x^3 + ax + b$ be an elliptic curve with non-zero discriminant, 
$\Delta(f)$, and real roots $e_1$, $e_2$, and $e_3$ of $f$,
with $e_1 < e_2 < e_3$.  Then we can
calculate a period basis, $\omega_1$ and $\omega_2$, as follows.

\begin{mydef}
Given $a,b \in \mathbb{R}$, let $a_0 = a$ and $b_0 = b$, and 
define the following recursion relation

\begin{align*}
a_{n+1} &= \frac{1}{2} \left( a_n + b_n \right),\\
b_{n+1} &= \sqrt{a_nb_n}.
\end{align*}
 
\end{mydef}

\begin{mythm}
In the previous definition, $a_n$ and $b_n$ converge to the same value.   We denote
this limit by $M(a,b) = a_{\infty} = b_{\infty}$ and call it the
\emph{arithmetic-geometric mean} of $a$ and $b$.
\end{mythm}

\begin{proof}
See \cite{Wash}.
\end{proof}

\begin{mythm}
Let $e_1 < e_2 < e_3$ be real numbers.  If we take $\Lambda =
\mathbb{Z}\omega_1 \oplus \mathbb{Z}\omega_2$ where,
\begin{align}
\omega_1 &= \frac{\pi i}{M\left(\sqrt{e_3-e_1}, \sqrt{e_2-e_1}\right)},\\
\omega_2 &= \frac{\pi}{M\left(\sqrt{e_3-e_1}, \sqrt{e_3-e_2}\right)},
\end{align}
then $E_\Lambda$ is given by $y^2 = (x-e_1)(x-e_2)(x-e_3)$.
\end{mythm}

\begin{proof}
See \cite{Wash}.
\end{proof}

\section{Group Properties}

The bijection between $(\mathbb{C}-\Lambda)/\Lambda$ and
solutions to $y^2 = x^3 + ax + b$ makes the points on an elliptic curve into a
group ``under addition.''  If we introduce a ``point at infinity,'' 
the identity element, $O$, and let $P = (x_P, y_P)$, $Q = (x_Q, y_Q)$ (where $P
\neq Q$ and $P, Q \neq O$), and define $P + Q = R$, where $R = (x_R, y_R)$, then
$s = \frac{y_P-y_Q}{x_P-x_Q}$, $x_R = s^2 - x_P - x_Q$, and $y_R = s(x_P-x_R) -
y_P$, (see \cite{Wikipedia}).
Geometrically, this is equivalent to intersecting the line through points
$P$ and $Q$ with the curve, and reflecting that intersection across the $x$-axis.

If $P = Q$, then the ``line through points $P$ and $Q$'' is the tangent line at
$P$.  If $P = O$, then $P + Q = Q$.

\begin{figure}[H]
\begin{center}
\includegraphics[scale=0.5]{EllipticCurve_800.png}
\caption{Point addition on an elliptic curve with a = -1, b = 0 \cite{Wolfram}.}
\end{center}
\end{figure}

\chapter{Division Polynomials}

If $n \in \mathbb{N}$ and $P$ is a point on an elliptic curve, define
\[ nP = \underbrace{P + \ldots + P}_{\text{n times}}. \]

We could perform this calculation the straight-forward way using the group law.  
This involves adding $P$ to itself $n$ times,
using the group addition we defined above.  This calculation can be sped up by adding $P$ to itself, then
adding that result to itself, and so on (repeated doubling).  So the first iteration would be $Q = P+P = 2P$.  The second 
iteration would be $R = Q + Q = 4P$.  So for example, $7P = R + Q + P$.  This methodology would be carried out until the desired result is 
achieved.  While this is a slight improvement, there is a quicker way if we
want to multiply many points by $n$.
This is where division polynomials comes into play.

\begin{mydef}
Given an elliptic curve, $y^2 = x^3 + ax + b$, let $f_n$ be defined as:
\begin{align*} 
f_1      &= 1 \\
f_2      &= 2 \\
f_3      &= 3x^4 + 6ax^2 + 12bx - a^2 \\
f_4      &= 4x^6 + 20ax^4 + 80bx^3 - 20a^2x^2 - 16abx - 32b^2 - 4a^3 \\
f_{2m}   &= f_m(f_{m+2}f_{m-1}^2 - f_{m-2}f_{m+1}^2)/2 \text{ for m} \geq 3\\
f_{4l+1} &= (x^3 + ax + b)^2 f_{2l+2}f_{2l} - f_{2l-1}f_{2l+1}^3 \text{ for l}\geq 1 \\
f_{4l+3} &= f_{2l+3}f_{2l+1}^3 -(x^3 + ax + b)^2f_{2l}f_{2l+2}^3\text{ for l}\geq 1
\end{align*}

Define the division polynomial, $\psi_n$ as:
\begin{align*}
\psi_n &= f_n \text{ if n is odd}, \\
\psi_n &= y f_n \text{ if n is even}
\end{align*}

See \cite{Div}.

\end{mydef}

\begin{mythm}
\[ nP = \left( \frac{\phi_n(x)}{\psi_n(x)^2}, \frac{\omega_n(x)}{\psi_n(x)^3}
\right) \]

where
 
\begin{align*} 
\phi_n &= x\psi_n^2 - \psi_{n+1}\psi_{n-1}, \\ 
\omega_n &= (4y)^{-1}(\psi_{n+2}\psi_{n-1}^2 - \psi_{n-2}\psi_{n+1}^2). 
\end{align*}
\end{mythm}

\begin{proof}
See \cite{Wash}.
\end{proof}

\begin{mythm}
If $P \neq O$, then $nP = O$ (point at infinity) if and only if $\psi_n(x) = 0$.
\end{mythm}

\begin{proof}
See \cite{Wash}.
\end{proof}

\chapter{Resultants}

Here's some background on resultants and discriminants.

\begin{mydef}
Let $f$ be a polynomial with roots $d_1, \ldots ,d_m$.  So we can write $f =
a_m(x-d_1)\dots(x-d_m)$ for some $a_m$.  The discriminant is defined as 
\[\Delta(f) = a_m^{2m-2}\myprod{1 \leq i,j \leq m}{i<j} \left(d_i -
d_j\right)^2\]
\end{mydef}

\begin{mydef}
Let $f = a_0 + a_1x + \ldots + a_mx^m$ and $g = b_0 + b_1x + \ldots + b_nx^n$ be two polynomials.  Then
the resultant of $f$ and $g$ with respect to $x$, denoted $res(f, g)$, is
defined as the determinant of the $(m+n) \times (m+n)$ Sylvester matrix $A$,
where \[ A = \left( \begin{array}{ccccccccc}
a_m    & a_{m-1}& a_{m-2}& \ldots & 0       & 0      & 0      \\ 
0      & a_m    & a_{m-1}& \ldots & 0       & 0      & 0      \\
\vdots & \vdots & \vdots & \ddots & \vdots  & \vdots & \vdots \\ 
0      & 0      & 0      & \ldots & a_1     & a_0    & 0      \\
0      & 0      & 0      & \ldots & a_2     & a_1    & a_0    \\
b_n    & b_{n-1}& b_{n-2}& \ldots & 0       & 0      & 0      \\ 
0      & b_n    & b_{n-1}& \ldots & 0       & 0      & 0      \\
\vdots & \vdots & \vdots & \ddots & \vdots  & \vdots & \vdots \\ 
0      & 0      & 0      & \ldots & b_1     & b_0    & 0      \\
0      & 0      & 0      & \ldots & b_2     & b_1    & b_0    \\
\end{array} \right). \]
\end{mydef}

\begin{mythm}
Let $f$ and $g$ be two polynomials with roots $d_1, \ldots ,d_m$ and $e_1,
\ldots ,e_n$, respectfully.  So we
can write $ f = a_m(x - d_1)\dots(x-d_m)$ and $g = b_n(x-e_1)\dots(x-e_n)$ for
some $a_m$, $b_n$.
Then
\[ res(f,g) = a_m^n b_n^m \myprod{1 \leq i \leq m}{1 \leq j \leq n} (d_i - e_j)
= a_m^n \prod_{1 \leq i \leq m}g(d_i) \]
\end{mythm}

\begin{proof}
See \cite{Alg}.
\end{proof}

Note that the resultant is zero if and only if $f$ and $g$ have a root in
common.  An equivalent definition of the discriminant is $\Delta(f) = res(f, f^\prime)$, 
which is zero if and only if $f$ has a multiple-root.

\chapter{Theta functions, Eta functions, Sigma functions}

\begin{mydef}
Given $a,b,z \in \mathbb{C}$, and $\tau \in \mathbb{H} = \left\{ x + iy | y > 0
\right\}$, we define the Riemann theta function with characteristic $\mych{a}{b}$ as

\begin{equation}
\mythv{a}{b}{z} = \sum_{n \in \mathbb{Z}}e^{i\pi(n+a)((n+a)\tau + 2(z+b)))}
\end{equation}
\end{mydef}

Four ``Jacobi,'' or ``elliptic,'' theta functions are often associated to $\tau$
and are given by
\begin{align}
  \vartheta_{00}(z;\tau) &= \mythv{0}{0}{z}, \\
  \vartheta_{01}(z;\tau) &= \mythv{0}{1/2}{z}, \\
  \vartheta_{10}(z;\tau) &= \mythv{1/2}{0}{z}, \\
  \vartheta_{11}(z;\tau) &= \mythv{1/2}{1/2}{z}.  
\end{align}

\begin{mydef}
We define the Dedekind eta function as
\[ \eta(\tau) = e^{\frac{\pi i \tau}{12}} \prod_{n=1}^{\infty} \left( 1-q^n
\right), \]

where $\tau \in \mathbb{H}$, and $q = e^{2 \pi i \tau}$.
\end{mydef}

\begin{mydef}
Let $\Lambda$ be a lattice in $\mathbb{C}$.  We define the sigma function
attached to $\Lambda$ as \[ \sigma(z, \Lambda) = z\myprod{\omega \in
\Lambda}{\omega \ne 0} \left( 1-\frac{z}{w} \right) e^{z/w + \frac{1}{2}
(z/w)^2}. \]
\end{mydef}

Now we list some identities related to the sigma, theta, and eta functions, some
of which I will verify later by computation, using explicit values.  See \cite{Acta}.

\begin{mythm}
The sigma, theta, and eta functions satisfy the following properties
\[
\eta(\tau)^{24} = \Delta(\tau) = g_2(\tau)^3 - 27g_3(\tau)^2,
\]
\[
\myprod{0 \le i,j < n}{(i,j) \ne (0,0)}\myth{\myvec{1/2 + a/n}{1/2 + b/n}}{0}
 = (-1)^{n-1} n \eta(\tau)^{n^2-1},
\]
\begin{align*}
\mythp{\myvec{1/2}{1/2}}{'}{0} &= -\pi \mythv{0}{0}{0} \mythv{1/2}{0}{0}
\mythv{0}{1/2}{0} \\ & = -2 \pi \eta(\tau)^3
\text{ (Jacobi's Derivative Formula)},
\end{align*}
\[
\mythv{a}{b}{0}^8 \myprod{0 \le i,j < 3}{(i,j) \ne (0,0)} 
\mythp{\myvec{a+i/3}{b+j/3}}{'}{0} = (-1)^{2a+1} 2^8 \frac{\pi^8}{3^5}
\eta(\tau)^{32},
\]
\[
\myprod{0 \le i,j < 4}{(i,j) \ne (0,0),(0,2),(2,0),(2,2)}
\mythp{\myvec{i/4}{j/4}}{'}{0} = \frac{-(\pi^{12})}{8} \eta(\tau)^{36}.
\]
\end{mythm}

\chapter{Resultants and Division Polynomials}
Now we're going to show how everything ties together.  From here on out, assume that an elliptic curve
is defined by $f = y^2 = x^3 + ax + b$, for some $a, b \in \mathbb{C}$.  Assume
this curve has a corresponding lattice $\Lambda$ with basis $\omega_1, \omega_2
\in \mathbb{C}$.  First note that the discriminant of $f$ is calculated to be
$-4a^3 - 27 b^2$, which must be non-zero.  That is,

 \[ \Delta = -16\left(4a^3 + 27b^2\right) \neq 0. \]

If we take the resultant of two division polynomials, we get some interesting
results.  The remainder of this chapter is a summary of a paper by Grant.  
See \cite{Grant} for details.

To begin, we need to define our derivation for the function field of $E$. 
Define $D = 2y \frac{d}{dx}$.  So, we can see that if $\psi$ is any polynomial
in $x$, and a prime denotes $\frac{d}{dx}$, then

\[ D^2\psi(x) = 4 \left( \psi''(x) y^2 + \psi'(x) \left( \frac{3x^2 + a}{2}
\right) \right). \]

Let $n$ be a positive odd integer.  First note that the $n^{th}$ division
polynomial may also be expressed by

\[ \psi_n(x) = \frac{\sigma(nz)}{\sigma(z)^{n^2}} \]
where $\psi_n$ is attached to $E_{\Lambda}$, $\Lambda = \mathbb{Z} \oplus
\mathbb{Z}\tau$, $x = \wp(z,\Lambda)$, $\tau \in \mathbb{H}$, and $\sigma(z) =
\sigma(z,\tau)$.

Also note that, by letting $\delta = \mych{1/2}{1/2}$,

\begin{equation}
\myth{\delta}{z} = ce^{dz^2}\sigma(z)
\end{equation}
where $c = \theta[\delta]'(0,\tau) = -2\pi\eta(\tau)^3$, and $d$ is independent
of $z$.

So we can substitute to get

\begin{equation}
\psi_n(x) = c^{n^2-1} \frac{\myth{\delta}{nz}}{\myth{\delta}{z}^{n^2}}
\end{equation}

Now let $a$, $b$ be integers, not both multiples of $n$.  Let $z_{a,b} =
\frac{a\tau+b}{n}$.  With our identifications, $D =
\frac{d}{dz}$, and applying $D^2$ to both sides of (6.2) gives

\begin{equation}
D^2\psi_n(x_{a,b}) =
c^{n^2-1}n^2\frac{\myth{\delta}{z_{a,b}}\mythp{\delta}{''}{nz_{a,b}} -
2n\mythp{\delta}{'}{z_{a,b}}}{\myth{\delta}{z_{a,b}}^{n^2+1}}
\end{equation}
where $x_{a,b} = \wp(z_{a,b})$ and the prime denotes $\frac{d}{dz}$.

This simplifies to

\[ D^2\psi_n(x_{a,b}) = 2c^{n^2}n^3(-1)^{(n-1)a+b+1}\frac{\mythp{\delta +
\myvec{a/n}{b/n}}{'}{0}} {\myth{\delta + \myvec{a/n}{b/n}}{0}^{n^2+1}}. \]

If we denote the resultant of $\psi_n(x))$ and $D^2\psi_n(x)$ as $R_n$, then by
Theorem 6,

\[ \frac{R_n^2}{n^{n^2+1}} = \myprod{0 \leq a,b < n}{(a,b) \neq
0}D^2\psi_n(x_{a,b}).
\]

Define

\[ Derivprod[n] = \myprod{0 \leq a,b < n}{(a,b) \neq
0}\mythp{\delta+\myvec{a/n}{b/n}}{'}{0}, \]

and

\[ Prod[n] = \myprod{0 \leq a,b < n}{(a,b) \neq
0}\myth{\delta+\myvec{a/n}{b/n}}{0}. \]

So we are left with

\[ \frac{R_n^2}{n^{n^2+1}} =
(2c^{n^2}n^3)^{n^2-1}(-1)^{\sum((n-1)a+b+1)}\frac{Derivprod[n]}{Prod[n]^{n^2+1}},
\]

which by Theorem 7 can be simplified to

\[ R_n^2 =
(-1)^{\frac{n-3}{2}}2^{n^4-1}\pi^{n^4-n^2}n^{3n^2-3}\eta(\tau)^{2n^4-3n^2+1}Derivprod[n]
\]
  
\chapter{Implementation}

At first, I began development using SAGE, an online math environment.  Although
basic implementations for elliptic curves, with division polynomials, were
provided, I quickly realized that I needed to implement my own functionality.

The language used for implementation was a combination of Python, SymPy, and
NumPy.  Python is a highly-extensible scripting language, with a wide variety of
applications.  It can be used as a
simple scripting language for quick results, or it can be used 
as a full-featured programming language, complete with classes, packages, etc.

NumPy is an extension to Python, which provides a variety of numerical
computation capabilities.  Although I use it for basic exponential calculations,
it does provide a wide variety of other numerical tools.

SymPy is an extension to Python, providing many features similar to
what Mathematica provides, include the ability to compute
expressions using symbolic constants.  This was useful for verification of
division polynomial calculations using indeterminate $a$'s and $b$'s (as in the
elliptic curve equation $y^2 = x^3 + ax + b$), instead of explicit values.  In
addition, SymPy provides the ability for arbitrary precision calculations.

The application used for development was Eclipse, which is a extremely
popular, highly extensible development environment.
Although traditionally used for Java development, Eclipse, along with the PyDev
plugin, is a powerful environment for Python development.  In addition, the
software was hosted remotely on Google's code repository.  The SVN plugin
permitted synchronization with the Google repository, allowing for convenient
distribution to any location.  To highlight
the variability of Eclipse, this document was written in LaTeX, using the Texclipse
plugin, and hosted on the Google repository as well.

\lstset{language=Python}

\section{EllipticCurve class}
The elliptic curve class is a basic software representation of a
Weierstrass elliptic curve, assuming the form $y^2 = x^3 + ax + b$.
The constructor takes 3 arguments, all of which are optional.  The basic form of
the constructor is EllipticCurve(r,s,isLattice).

\begin{lstlisting}
class EllipticCurve:
    def __init__(self, r=None, s=None, isLattice=False):
\end{lstlisting}

Essentially there are three use cases. 

If either $r$ or $s$ is
undefined, isLattice is implicitly false, and the elliptic curve is initialized
with $a$ and $b$ as indeterminates.
The elliptic curve $y^2 = x^3 + ax + b$ is created with symbolic constants
$a$ and $b$ (using SymPy), and the division polynomials returned are symbolic,
with indeterminate $a$'s and $b$'s.

\begin{lstlisting}
        if r == None or s == None:
            self.a = sympy.abc.a
            self.b = sympy.abc.b
            return
\end{lstlisting}

If isLattice is true, and both $r$ and $s$ are provided, then $r$ and $s$
represent a basis for a lattice, $\omega_1 = r$, $\omega_2 = s$, and $\Lambda$
is the $\mathbb{Z}$-span of $\omega_1$ and $\omega_2$.  The values $g_2$ and
$g_3$ are calculated using (see \cite{Handbook} and \cite{Wikipedia})

\begin{align} 
  g_2(\omega_1,\omega_2) &= \frac{2}{3} \left(\frac{\pi}{\omega_1}\right)^4
  \left(\vartheta_{10}(0,\tau)^8 + \vartheta_{00}(0,\tau)^8 +
  \vartheta_{01}(0,\tau)^8\right) \\
   g_3(\omega_1,\omega_2) &= \left(\frac{\pi}{\omega_1}\right)^6 \left[ 
  \frac{8}{27} \left( \vartheta_{10}(0,\tau)^{12} + \vartheta_{00}(0,\tau)^{12}
  \right) \right. - \\
  &\left. \frac{4}{9}\left(  \vartheta_{10}(0,\tau)^4 + \vartheta_{00}(0,\tau)^4
  \right) \vartheta_{10}(0,\tau)^4 \vartheta_{00}(0,\tau)^4 \right]
\end{align}

Finally $a$ and $b$ in $y^2 = x^3 + ax + b$ is calculated using

\begin{align}
  a &= -\frac{g_2}{4} \\
  b &= -\frac{g_3}{4}
\end{align}

\begin{lstlisting}
        if isLattice == True:
            self.w1 = r
            self.w2 = s
            if ((self.w2/self.w1).imag < 0.0):
                self.w1 = r
                self.w2 = s

            self.tau = self.w2/self.w1
            
            g2 = (2.0/3.0) * (pi/self.w1)**4 * (jacobiTheta(1,0,0.0,self.tau)**8 + jacobiTheta(0,0,0.0,self.tau)**8 + jacobiTheta(0,1,0.0,self.tau)**8)
            g3 = (pi/self.w1)**6 * ((8.0/27.0) * (jacobiTheta(1,0,0.0,self.tau)**12 + jacobiTheta(0,0,0.0,self.tau)**12) - (4.0/9.0)*(jacobiTheta(1,0,0.0,self.tau)**4 + jacobiTheta(0,0,0.0,self.tau)**4) * jacobiTheta(1,0,0.0,self.tau)**4 * jacobiTheta(0,0,0.0,self.tau)**4)
            self.a = (-1.0/4.0) * g2
            self.b = (-1.0/4.0) * g3
\end{lstlisting}

If $r$ and $s$ are both defined, and isLattice is false (the default value),
then the elliptic curve $y^2 = x^3 + ax + b$ is initialized with $a = r$ and $b
= s$.
In addition, the roots of the polynomial are calculated using NumPy's polynomial
root calculator.  Finally, the periods (the basis of the associated lattice) are
calculated using the arithmetic-geometric mean.  Note that because the
arithmetic-geometric mean, as implemented, works only on real values, we must
assure that the roots of the polynomial are real.  We cannot construct an arbitrary elliptic curve, with any $a,b \in \mathbb{C}$.  

\begin{lstlisting}
        else:
            self.a = r
            self.b = s
            roots = sorted(poly1d([1,0,a,b]).r)
            self.e1 = roots[0]
            self.e2 = roots[1]
            self.e3 = roots[2]
             
            self.w1 = pi * (0.0+1.0j) / \
                self.arithGeomMean(sqrt(self.e3-self.e1),\
                                   sqrt(self.e2-self.e1))
            self.w2 = pi / \
                self.arithGeomMean(sqrt(self.e3-self.e1),\
                                   sqrt(self.e3-self.e2))
            self.tau = self.w1/self.w2
            if (self.tau.imag < 0.0):
                self.tau = 1/self.tau
\end{lstlisting}

The arithmetic geometric mean is defined by the recursion relations 
\begin{align}
a_{n+1} &= \frac{1}{2} \left( a_n + b_n \right)\\
b_{n+1} &= \sqrt{a_nb_n}
\end{align}

As noted above, this implementation works only for real values.

\begin{lstlisting}
    def arithGeomMean(self,a,b):
        anew = a
        bnew = b
        tol = 0.0000000000001
        while(anew - bnew >= tol):
            aprev = anew
            bprev = bnew
            anew = 0.5 * (aprev+bprev)
            bnew = sqrt(aprev*bprev)
        return anew
\end{lstlisting}

The following are the division polynomial implementations.

\begin{lstlisting}
    def divPoly(self, n):
        if n%2 == 0:
            return self.divPolyHelp(n)*Poly([1,0],y)
        else:
            return self.divPolyHelp(n)
        
    def divPolyHelp(self,n):
        if n == 1:
            return(Poly([1], x))
        elif n == 2:
            return(Poly([2], x))
        elif n == 3:
            return(Poly([3, 0, 6*self.a, 12*self.b, -self.a**2], x))
        elif n == 4:
            return(Poly([4,0, 20*self.a, 80*self.b, \
                               -20*self.a**2, -16*self.a*self.b, \
                               -32*self.b**2 - 4*self.a**3], x))
        elif n%2 == 0:
            m = n/2
            return(Poly(self.divPolyHelp(m) * (self.divPolyHelp(m+2)*\
                        self.divPolyHelp(m-1)**2 - self.divPolyHelp(m-2)*\
                        self.divPolyHelp(m+1)**2)/2))
        elif n%4 == 1:
            l = n/4
            return(Poly([1,0,self.a,self.b], x)**2 * \
                   self.divPolyHelp(2*l + 2) * self.divPolyHelp(2*l)**3 - \
                   self.divPolyHelp(2*l-1)*self.divPolyHelp(2*l+1)**3)
        elif n%4 == 3:
            l = n/4
            return(self.divPolyHelp(2*l+3)*self.divPolyHelp(2*l+1)**3 - \
                   Poly([1,0,self.a,self.b], x)**2 * \
                   self.divPolyHelp(2*l)*self.divPolyHelp(2*l+2)**3)
        
\end{lstlisting}

For the second derivative, we define 
\[ D^2\psi(x) = 4 \left( \psi''(x) y^2 + \psi'(x) \left( \frac{3x^2 + a}{2} \right) \right) \]

\begin{lstlisting}
    def secondDeriv(self, f):
        return 4*(f.diff(x).diff(x) * \
                  Poly([1, 0, self.a, self.b], x) \
                  + f.diff(x) * Poly([3, 0, self.a], x)/2.0)
\end{lstlisting}

The Weierstrass $\wp$-function is defined as
\[ \wp(z) = \frac{1}{z^2} + \sum_{\substack{\omega \in \Lambda\\\omega \ne 0}} 
\left( \frac{1}{(z-\omega)^2} - \frac{1}{\omega^2} \right) \]
\begin{lstlisting}
    def pFunc(self, z):
        z = z/self.w1
        val = pi**2 * jacobiTheta(0,0,0.0,self.tau)**2 * \
            jacobiTheta(1,0,0.0,self.tau)**2 * \
            jacobiTheta(0,1,z,self.tau)**2 / \
            jacobiTheta(1,1,z,self.tau)**2 - pi**2/3.0 * \
            (jacobiTheta(0,0,0.0,self.tau)**4 + \
             jacobiTheta(1,0,0.0,self.tau)**4)
        return val/(self.w1**2)
\end{lstlisting}

Return the discriminant, 

\[ \Delta = -16\left(4a^3 + 27b^2\right) \]

\begin{lstlisting}                    
    def discrim(self):
        return -16*(4*self.a**3 + 27*self.b**2)
\end{lstlisting}

\section{ThetaFunction module}
This file provides functionality related to theta functions with
characteristic.
This includes the theta function, and its first and second derivatives.  A Jacobi theta function
is provided for the $\vartheta_{00}$, $\vartheta_{01}$, $\vartheta_{10}$, and
$\vartheta_{11}$ forms of the theta function.  The mappings between the Jacobi
theta and the theta function with characteristic is provided earlier in this
paper.

\begin{lstlisting}
def theta(a,b,z,t):
    temp = (0 + 0j)
    temp = temp + helperFunc(a,b,z,t,0)
    for n in range(1,1000):
        temp = temp + helperFunc(a,b,z,t,n)
        temp = temp + helperFunc(a,b,z,t,-n)
    return temp

def jacobiTheta(n,m,z,t):
    a = 0.5 if n == 1 else 0.0
    b = 0.5 if m == 1 else 0.0
    return theta(a,b,z,t)

def thetaDeriv(a,b,z,t):
    temp = (0 + 0j)
    temp = temp + 2*pi*(0+1j)*a*helperFunc(a,b,z,t,0)
    for n in range(1,1000):
        temp = temp + 2*pi*(0+1j)*(n+a)*helperFunc(a,b,z,t,n)
        temp = temp + 2*pi*(0+1j)*(-n+a)*helperFunc(a,b,z,t,-n)
    return temp

def thetaSecondDeriv(a,b,z,t):
    temp = (0 + 0j)
    temp = temp + ((2*pi*(0+1j)*a)**2)*helperFunc(a,b,z,t,0)
    for n in range(1,1000):
        temp = temp + ((2*pi*(0+1j)*(n+a))**2)*helperFunc(a,b,z,t,n)
        temp = temp + ((2*pi*(0+1j)*(-n+a))**2)*helperFunc(a,b,z,t,-n)
    return temp

def helperFunc(a,b,z,t,n):
    return exp(pi * (0 + 1j) * (n+a) * ((n+a)*t + 2*(z+b)))
\end{lstlisting}

\section{SigmaFunction module}
This file provides the sigma function,
\[ \sigma(z, \Lambda) = z\myprod{\omega \in \Lambda}{\omega \ne (0,0)} \left(
1-\frac{z}{w} \right) e^{z/w + \frac{1}{2} (z/w)^2} \]

\begin{lstlisting}
def sigma(z,w1,w2):
    res = z
    for m in range(-500,501):
        for n in range(-500,501):
            if m == 0 and n == 0:
                continue
            w = m*w1 + n*w2
            res *= ((1-z/w)*exp(z/w + 0.5*(z/w)**2))
    return res
\end{lstlisting}

\section{EtaFunction module}
This file provides the eta function,

\[ \eta(\tau) = e^{\frac{\pi i \tau}{12}} \prod_{n=1}^{\infty} \left( 1-q^n \right) \]

\begin{lstlisting}
def eta(t):
    q = exp(2*pi*(0.0+1.0j)*t)
    res = exp((0.0+1.0j)*pi*t/12.0)
    for n in range(1000):
        res = res * (1-power(q,n+1))
    return res
\end{lstlisting}

\section{ThetaProds module}
This file provides the derivProd and prod functions, as used in \cite{Acta} and \cite{Grant}.  The derivProd and prod functions are defined as:

\begin{align} 
Derivprod[n] &= \myprod{a,b < n}{(a,b) \neq
0}\mythp{\delta+\myvec{a/n}{b/n}}{'}{0}\\
Prod[n]      &= \myprod{a,b < n}{(a,b) \neq 0}\myth{\delta+\myvec{a/n}{b/n}}{0}
\end{align}

\begin{lstlisting}
def derivProd(n, t, a = 0.5, b = 0.5, excl = [(0,0)]):
    res = 1.0
    for i in range(n):
        for j in range(n):
            if (i,j) in excl:
                continue
            res = res * thetaDeriv(a + float(i)/n, b + float(j)/n, 0.0, t)
    return res
 
def prod(n,t):
   res = 1.0
   for i in range(n):
       for j in range(n):
           if i == 0 and j == 0:
               continue
           res = res * theta(0.5 + float(i)/n, 0.5 + float(j)/n, 0.0, t)
   return res
\end{lstlisting}

As you can see, the derivProd function takes five parameters.  The $a$, $b$, and $excl$ parameters have default values of $0.5$, $0.5$, and $[(0,0)]$, 
respectfully.  The $excl$ parameter defines what $(a,b)$ value pairs to not use in the calculation.  For example, the derivProd4 function as
defined in \cite{Acta} is implemented as $derivProd(4, e.tau, 0.0, 0.0, [(2,0),(0,2),(0,0),(2,2)])$.



\chapter{Computational Results}

\newcommand{\shellcmd}[1]{\noindent\indent\indent\texttt{\footnotesize>>
#1}\\}
\newcommand{\rescmd}[1]{\noindent\indent\indent\texttt{\footnotesize#1}\\}
\newcommand{\dotscmd}[1]{\noindent\indent\indent\texttt{\footnotesize...
#1}\\}

For the following results, we will attempt to initialize the elliptic curves
with a wide variety of parameters, to make a good faith effort of proving
correctness.
Keep in mind that the elliptic curve could not be initialized from the $a$'s and $b$'s
(non-lattice), where the elliptic curve would result in imaginary roots
(the arithmetic-geometric mean routine only works with real values).  Also,
for calculations that involve division polynomials in relation to theta
functions (i.e. the right hand side depended only on $\tau$, instead of
$\omega_1$ and $\omega_2$), we were limited to elliptic curves with a lattice
basis of the form $\omega_1$, $\omega_2$, where $\omega_1 = 1$.  When
referenced, $\delta = \left[ \myvec{1/2}{1/2} \right] $.

For the first test, we will attempt to check that 

\[ Prod[n] = \myprod{0 \leq a,b < n}{(a,b) \neq
0}\myth{\delta+\myvec{a/n}{b/n}}{0} = (-1)^{n-1} n \eta(\tau)^{n^2-1}\]

\shellcmd{e = EllipticCurve(1,2) }
\shellcmd{n = 7}
\shellcmd{print prod(n, e.tau)}
\rescmd{(-1.0205452215e-05+4.75920707408e-06j)}
\shellcmd{print (-1)**(n-1) * n * eta(e.tau)**(n**2-1)}
\rescmd{(-1.0205452215e-05+4.75920707408e-06j)}

Next we will check that 
\[ \mythp{\delta}{'}{0} = -\pi \mythv{0}{0}{0} \mythv{1/2}{0}{0}
\mythv{0}{1/2}{0} = -2 \pi \eta(\tau)^3 \]

This is also known as Jacobi's derivative formula.
\\

\shellcmd{e = EllipticCurve((0.0+1.0j),1.0, True)}
\shellcmd{print thetaDeriv(0.5,0.5,0.0,e.tau)}
\rescmd{(-2.84869460399+1.82761721602e-33j)}
\shellcmd{print -pi * theta(0.0,0.0,0.0,e.tau) * theta(0.5,0.0,0.0,e.tau) * \\
theta(0.0,0.5,0.0,e.tau)}
\rescmd{(-2.84869460399+0j)}
\shellcmd{print -2 * pi * eta(e.tau)**3}
\rescmd{(-2.84869460399+0j)}

Now we will check that 
\begin{align*}
\mythv{a}{b}{0}^8 \myprod{0 \le i,j < 3}{(i,j) \ne (0,0)} 
\mythp{\myvec{a+i/3}{b+j/3}}{'}{0} &= (-1)^{2a+1} 2^8 \frac{\pi^8}{3^5}
\eta(\tau)^{32} \\
\myprod{0 \le i,j < 4}{(i,j) \ne (0,0),(0,2),(2,0),(2,2)}
\mythp{\myvec{i/4}{j/4}}{'}{0} &= \frac{-(\pi^{12})}{8} \eta(\tau)^{36}
\end{align*}

Where $(a,b) = (0,0)$, $(0,1/2)$, or $(1/2,0)$.
\\

\shellcmd{e = EllipticCurve(-4,5)}
\shellcmd{for (a,b) in [(0,0),(0,0.5),(0.5,0)]:}
\dotscmd{    print a,b}
\dotscmd{    print theta(a,b,0.0,e.tau)**8 * derivProd(3,e.tau,a,b)}
\dotscmd{    print ((-1)**(2*a+1) * 2**8 * pi**8/(3.0**5)) * eta(e.tau)**32}
\dotscmd{}
\rescmd{0 0}
\rescmd{(-2.16506165686+1.29297718511e-15j)}
\rescmd{(-2.16506165686+0j)}
\rescmd{0 0.5}
\rescmd{(-2.16506165686+2.61289139492e-15j)}
\rescmd{(-2.16506165686+0j)}
\rescmd{0.5 0}
\rescmd{(2.16506165686-2.9091986665e-15j)}
\rescmd{(2.16506165686+0j)}
\shellcmd{print derivProd(4, e.tau, 0.0, 0.0, [(2,0),(0,2),(0,0),(2,2)])}
\rescmd{(-8.71567004001+8.881784197e-16j)}
\shellcmd{print (-(pi**12)/8.0) * eta(e.tau)**36}
\rescmd{(-8.71567004001+0j)}

This test will verify that 
\[ \psi_n = \frac{\sigma(nz)}{\sigma(z)^{n^2}} =
c^{n^2-1}\frac{\myth{\delta}{nz}}{\myth{\delta}{z}^{n^2}} \text{ , where } c =
-2\pi\eta(\tau)^3 \].
 
\shellcmd{e = EllipticCurve((0.5+1.25j),1.0,True)}
\shellcmd{z = (1.25+0.75j)}
\shellcmd{n = 7}
\shellcmd{f = e.divPoly(n)}
\shellcmd{print sigma(n*z, e.w1, e.w2)/(sigma(z, e.w1, e.w2)**(n**2))}
\rescmd{(-1.2638839763e+18-3.16870903571e+18j)}
\shellcmd{print f.eval(e.pFunc2(z))}
\rescmd{-1.2664020260877e+18 - 3.17530702355144e+18*I}
\shellcmd{c = -2 * pi * eta(e.tau)**3}
\shellcmd{print c**(n**2-1) * theta(0.5,0.5,n*z,e.tau)/\\(theta(0.5,0.5,z,e.tau)**(n**2))}
\rescmd{(-1.26639878842e+18-3.17522471162e+18j)}

Now we verify
\[
D^2\psi_n(x_{a,b}) =
c^{n^2-1}n^2\frac{\myth{\delta}{z_{a,b}}\mythp{\delta}{''}{nz_{a,b}} -
2n\mythp{\delta}{'}{z_{a,b}}}{\myth{\delta}{z_{a,b}}^{n^2+1}}
\]

\shellcmd{e = EllipticCurve((0.25+1.75j),1,True)}
\shellcmd{n = 5}
\shellcmd{a = 6}
\shellcmd{b = 7}
\shellcmd{f = e.divPoly(n)}
\shellcmd{g = e.secondDeriv(f)}
\shellcmd{c = -2 * pi * eta(e.tau)**3}
\shellcmd{z = (a*e.tau+b)/float(n)}
\shellcmd{fac1 = n**2 * c**(n**2-1)}
\shellcmd{tz = theta(0.5,0.5,z,e.tau)}
\shellcmd{tpz = thetaDeriv(0.5,0.5,z,e.tau)}
\shellcmd{tpnz = thetaDeriv(0.5,0.5,n*z,e.tau)}
\shellcmd{tppnz = thetaSecondDeriv(0.5,0.5,n*z,e.tau)}
\shellcmd{print g.eval(e.pFunc2(z))}
\rescmd{406358220709.026 + 520207245995.588*I}
\shellcmd{print fac1 * (- 2*n*tpz*tpnz + tppnz*tz)/(tz**(n**2+1))}
\rescmd{(406346558430+520203366532j)}

Now we verify

\[ D^2\psi_n(x_{a,b}) = 2c^{n^2}n^3(-1)^{(n-1)a+b+1}\frac{\mythp{\delta +
\myvec{a/n}{b/n}}{'}{0}} {\myth{\delta + \myvec{a/n}{b/n}}{0}^{n^2+1}} \]

\shellcmd{e = EllipticCurve((5.5+0.2j),1,True)}
\shellcmd{f = e.divPoly(n)}
\shellcmd{g = e.secondDeriv(f)}
\shellcmd{z = (a * e.tau + b)/float(n)}
\shellcmd{c = -2 * pi * eta(e.tau)**3}
\shellcmd{print 2 * c**(n**2) * n**3 * (-1)**(b+1) * thetaDeriv(0.5 + float(a)/n, \\0.5 + float(b)/n, 0.0, e.tau)/(theta(0.5 + float(a)/n, 0.5 + float(b)/n, 0.0, \\e.tau)**(n**2+1))}
\rescmd{(5.50700309786e+25-5.81172844367e+12j)}
\shellcmd{print complex(g.eval(e.pFunc2(z)))}
\rescmd{(5.5070054085e+25-2.15517056473e+18j)}

This final test of Grant's paper will verify that  

\[ 
\frac{R_n^2}{n^{n^2+1}} = (-1)^{\frac{n-3}{2}}2^{n^4-1}\pi^{n^4-n^2}n^{2n^2-4}\eta(\tau)^{2n^4-3n^2+1}Derivprod[n]
\]

\shellcmd{e = EllipticCurve((0.25+0.75j),1,True)}
\shellcmd{n = 7}
\shellcmd{f = e.divPoly(n)}
\shellcmd{g = e.secondDeriv(f)}
\shellcmd{print ((-1)**(n**2 * (n-1)/2 - 1) * Int(2)**(n**4-1) *\\
(Int(1)*pi)**(n**4-n**2)*(Int(1)*eta(e.tau))**(2*n**4-3*n**2+1)*\\
n**(2*n**2-4)*derivProd(n,e.tau)).evalf(10)}
\rescmd{-2.719526684e+1579 - 3.790087443e+1578*I}
\shellcmd{print ((f.resultant(g)**2)/(n**(n**2+1))).evalf(10)}
\rescmd{-2.719529682e+1579 - 3.78973944e+1578*I}

As we saw before, a division polynomial can be defined as $\psi_n = f_n$ if $n$
is odd, or $\psi_n = yf_n$ if $n$ is even.  This test will verify that
$res(\frac{f_n}{2},f_m)$ is a power of the discriminant (when $n$ is odd, $m$
is even) and $res(f_n, f_m)$ is a power of the discriminant (when $n$ and $m$
are odd).
The first case verifies a scenario when $n$ and $m$ are odd.

\shellcmd{e = EllipticCurve(-25,10)}
\shellcmd{fn = e.divPoly(9)}
\shellcmd{fm = e.divPoly(11)}
\shellcmd{res = fn.resultant(fm)}
\shellcmd{print res}
\rescmd{213042700113100244253453205893975276534812791049550171975960232262256435185\\
		3341589915627702667554435092936784813615943684831477708573860889361842782631775926\\
		2879279394396627174700966430396128839577698446275854008019673780263618071129929432\\
		6216438865223702361632791602140754277417262920611958705224773374644525426299965555\\
		8878803386403603207543763188902743141206220072597125050626751718339644497738834331\\
		7194730188462198839115239654821201401552948136233352939003939981701894179751624492\\
		9625063848574835556451033471657699869928953885207190625721586301691221636593279455\\
		0902379786661632617304316144742326631862491461799951703098402920201405786475390243\\
		7022310761550132792495922370100613480083674046486673490027668436168216600380771489\\
		0164781902591317874537713174280413638639897119384736442748417861615238101932431002\\
		2030631369486936452250631767718944088406553818855480161621023972904587472526217885\\
		4819067494106370717876391456857589191463925620699266093753005072309356089623102360\\
		9640471868351282147887306644115541865464648610843938164949135265492894988510639218\\
		0082906734544545934352447122630060939901775971644364649973673977542117851913140625\\
		6844140811388706794521614311658773504354300521615346875416424059568462091343161893\\
		2924368557737803656142650004976200015381951180902860021251724340898408909766756082\\
		2182128630025181332809896430427435646808498335455798081785514491740255143621391222\\
		0894965703164017937069145722225800243652635477708323414612530188121812041549947563\\
		7576752837125762530655356615054163176824641709057652603993621146733005005355105942\\
		4863848900034383732139291006893924559093760000000000000000000000000000000000000000\\
		0000000000000000000000000000000000000000000000000000000000000000000000000000000000\\
		0000000000000000000000000000000000000000000000000000000000000000000000000000000000\\
		0000000000000000000000000000000000000000000000000000000000000000000000000000000000\\
		0000000000000000000000000000000000000000000000000000000000000000000000000000000000\\
		0000000000000000000000000000000000000000000000000000000000000000000000000000000000\\
		0000000000000000000000000000000000000000000000000000000000000000000000000000000000\\
		0000000000000000000000000000000000000000000000000000000000000000000000000000000000\\
		0000000000000000000000000000000000000000000000000000000000000000000000000000000000\\
		0000000000000000000000000000000000000000000000000000000000000000000000000000000000\\
		0000000000000000000000}
\shellcmd{d = e.discrim()}
\shellcmd{print d}
\rescmd{956800}
\shellcmd{print res/(d**400)}
\rescmd{1}

This second case verifies a scenario when $n$ is even, $m$ is odd.

\shellcmd{e = EllipticCurve(-10,1)}
\shellcmd{fn = e.divPoly(8)}
\shellcmd{print fn}
\rescmd{Poly(8*x**30*y - 5840*x**28*y + 10016*x**27*y - 1087200*x**26*y + \\
187200*x**25*y + 96423552*x**24*y - 63648000*x**23*y - 2154560640*x**22*y - \\
536520192*x**21*y + 114237484800*x**20*y - 33914764800*x**19*y - \\
2408373805568*x**18*y + 1849041331200*x**17*y + 12245121085440*x**16*y - \\
20164971528192*x**15*y + 131897835264000*x**14*y + 69685511454720*x**13*y - \\
2441452139737088*x**12*y + 951883928985600*x**11*y + 11183441076510720*x**10*y - \\
6034621166927872*x**9*y - 19575319341465600*x**8*y + 314866672926720*x**7*y + \\
96414990066155520*x**6*y - 31341797169561600*x**5*y - 103370376850636800*x**4*y + \\
48701854875779072*x**3*y - 64064105254092800*x**2*y + 9483089722736640*x*y + \\
7301184751665152*y, x, y, domain='ZZ')}
\shellcmd{(fn, r) = div(fn, 2*y) fm = e.divPoly(7)}
\shellcmd{res = fn.resultant(fm)}
\shellcmd{d = e.discrim()}
\shellcmd{res.eval(1)}
\rescmd{244695355025057726866606179876151727893005016452444631462460853698539506247\\
		5760311422264929747731352258241802850699658489528586554902831705922689478206260219\\
		0390963439361155731116315300173309254437618565099311773125564668172278642893822942\\
		8510383842432787792208572102783107072618473683470584811842032964145956415045193009\\
		7313247670428375126103545123450417122757739681142330369831160154514430221314114364\\
		0301855049295051569952892934462909231463997437127199224535200880412137491517058523\\
		0861130169237342749875574541529878549845861098613192747393788592371215957811160868\\
		2107109376}
\shellcmd{print res/(d**120)}
\rescmd{1}

Finally, we show a result computing a resultant of two division polynomials,
with indeterminate $a$'s and $b$'s. 
\\

\shellcmd{e = EllipticCurve()}
\shellcmd{fn = e.divPoly(3)}
\shellcmd{print fn}
\rescmd{Poly(3*x**4 + 6*a*x**2 + 12*b*x - a**2, x, domain='ZZ[a,b]')}
\shellcmd{fm = e.divPoly(7)}
\shellcmd{print fm}
\rescmd{Poly(7*x**24 + 308*a*x**22 + 3944*b*x**21 - 2954*a**2*x**20 - 112*a*b*x**19 \\+
(-19852*a**3 - 42896*b**2)*x**18 - 92568*a**2*b*x**17 + \\
(-35231*a**4 - 571872*a*b**2)*x**16 + (-31808*a**3*b - 829696*b**3)*x**15 + \\
(-82264*a**5 - 615360*a**2*b**2)*x**14 + (-161840*a**4*b - 2132480*a*b**3)*x**13 + \\
(-111916*a**6 - 297472*a**3*b**2 - 928256*b**4)*x**12 + \\
(-608160*a**5*b - 2603776*a**2*b**3)*x**11 + \\
(-42168*a**7 - 1192800*a**4*b**2 - 3293696*a*b**4)*x**10 + \\
(-425712*a**6*b - 3727360*a**3*b**3 - 1555456*b**5)*x**9 + \\
(15673*a**8 - 831936*a**5*b**2 - 7069440*a**2*b**4)*x**8 + \\
(-53824*a**7*b - 1314560*a**4*b**3 - 7127040*a*b**5)*x**7 + \\
(14756*a**9 - 190400*a**6*b**2 - 2293760*a**3*b**4 - 2809856*b**6)*x**6 + \\
(57288*a**8*b - 168448*a**5*b**3 - 3698688*a**2*b**5)*x**5 + \\
(1302*a**10 + 134400*a**7*b**2 + 394240*a**4*b**4 - 3039232*a*b**6)*x**4 + \\
(1680*a**9*b + 152320*a**6*b**3 + 831488*a**3*b**5 - 802816*b**7)*x**3 + \\
(196*a**11 + 3696*a**8*b**2 + 96768*a**5*b**4 + 544768*a**2*b**6)*x**2 + \\
(392*a**10*b + 7168*a**7*b**3 + 64512*a**4*b**5 + 229376*a*b**7)*x - a**12 + \\
160*a**9*b**2 + 3328*a**6*b**4 + 24576*a**3*b**6 + 65536*b**8, x, domain='ZZ[a,b]')}
\shellcmd{res = fn.resultant(fm)}
\shellcmd{print res}
\rescmd{79228162514264337593543950336*a**48 + \\
8556641551540548460102746636288*a**45*b**2 + \\
433179978546740265792701548462080*a**42*b**4 + \\
13645169324222318372470098776555520*a**39*b**6 + \\
299340902050127109296062791910686720*a**36*b**8 + \\
4849322613212059170596217228953124864*a**33*b**10 + \\
60010367338499232236128188208294920192*a**30*b**12 + \\ 
578671399335528310848378957722843873280*a**27*b**14 + \\
4394285938704168110504877710207845662720*a**24*b**16 + \\
26365715632225008663029266261247073976320*a**21*b**18 + \\
124578006362263165932813283084392424538112*a**18*b**20 + \\
458673568879241656388994360447081199435776*a**15*b**22 + \\
1290019412472867158594046638757415873413120*a**12*b**24 + \\
2679271087443647175541481480496171429396480*a**9*b**26 + \\
3875374251480989664622499998574819388948480*a**6*b**28 + \\
3487836826332890698160249998717337450053632*a**3*b**30 + \\
1471431161109188263286355468208876736741376*b**32}
\shellcmd{res.factor()}
\rescmd{18446744073709551616*(4*a**3 + 27*b**2)**16}
\shellcmd{d = e.discrim()}
\shellcmd{print d}
\rescmd{-64*a**3 - 432*b**2}
\shellcmd{16**16}
\rescmd{18446744073709551616L}

\newpage

\bibliographystyle{alpha} 
\bibliography{ellipticcurve}

\end{document}